{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 数据准备"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "相关系数矩阵为：\n",
      "        x1    x2    x3    x4    x5    x6    x7    x8    x9   x10   x11   x12  \\\n",
      "x1   1.00  0.95  0.95  0.97  0.97  0.99  0.95  0.97  0.98  0.98 -0.29  0.94   \n",
      "x2   0.95  1.00  1.00  0.99  0.99  0.92  0.99  0.99  0.98  0.98 -0.13  0.89   \n",
      "x3   0.95  1.00  1.00  0.99  0.99  0.92  1.00  0.99  0.98  0.99 -0.15  0.89   \n",
      "x4   0.97  0.99  0.99  1.00  1.00  0.95  0.99  1.00  0.99  1.00 -0.19  0.91   \n",
      "x5   0.97  0.99  0.99  1.00  1.00  0.95  0.99  1.00  0.99  1.00 -0.18  0.90   \n",
      "x6   0.99  0.92  0.92  0.95  0.95  1.00  0.93  0.95  0.97  0.96 -0.34  0.95   \n",
      "x7   0.95  0.99  1.00  0.99  0.99  0.93  1.00  0.99  0.98  0.99 -0.15  0.89   \n",
      "x8   0.97  0.99  0.99  1.00  1.00  0.95  0.99  1.00  0.99  1.00 -0.15  0.90   \n",
      "x9   0.98  0.98  0.98  0.99  0.99  0.97  0.98  0.99  1.00  0.99 -0.23  0.91   \n",
      "x10  0.98  0.98  0.99  1.00  1.00  0.96  0.99  1.00  0.99  1.00 -0.17  0.90   \n",
      "x11 -0.29 -0.13 -0.15 -0.19 -0.18 -0.34 -0.15 -0.15 -0.23 -0.17  1.00 -0.43   \n",
      "x12  0.94  0.89  0.89  0.91  0.90  0.95  0.89  0.90  0.91  0.90 -0.43  1.00   \n",
      "x13  0.96  1.00  1.00  1.00  0.99  0.94  1.00  1.00  0.99  0.99 -0.16  0.90   \n",
      "y    0.94  0.98  0.99  0.99  0.99  0.91  0.99  0.99  0.98  0.99 -0.12  0.87   \n",
      "\n",
      "      x13     y  \n",
      "x1   0.96  0.94  \n",
      "x2   1.00  0.98  \n",
      "x3   1.00  0.99  \n",
      "x4   1.00  0.99  \n",
      "x5   0.99  0.99  \n",
      "x6   0.94  0.91  \n",
      "x7   1.00  0.99  \n",
      "x8   1.00  0.99  \n",
      "x9   0.99  0.98  \n",
      "x10  0.99  0.99  \n",
      "x11 -0.16 -0.12  \n",
      "x12  0.90  0.87  \n",
      "x13  1.00  0.99  \n",
      "y    0.99  1.00  \n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "data = pd.read_csv('../data/data.csv')  # 读取数据\n",
    "# 保留两位小数，并将结果保存为’.csv’文件\n",
    "np.round(data.corr(method = 'pearson'), 2).to_csv('../tmp/data_cor.csv')\n",
    "print('相关系数矩阵为：\\n', np.round(data.corr(method = 'pearson'), 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 特征工程"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "相关系数为： [ -1.80000000e-04  -0.00000000e+00   1.24140000e-01  -1.03100000e-02\n",
      "   6.54000000e-02   1.20000000e-04   3.17410000e-01   3.49000000e-02\n",
      "  -0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00\n",
      "  -4.03000000e-02]\n",
      "相关系数非零个数为： 8\n",
      "相关系数是否为零： [ True False  True  True  True  True  True  True False False False False\n",
      "  True]\n",
      "输出数据的维度为： (20, 8)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "D:\\DELL\\AppData\\Anaconda3\\lib\\site-packages\\sklearn\\linear_model\\coordinate_descent.py:491: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations. Fitting data with very small alpha may cause precision problems.\n",
      "  ConvergenceWarning)\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from sklearn.linear_model import Lasso\n",
    "\n",
    "data = pd.read_csv('../data/data.csv')  # 读取数据\n",
    "# 调用Lasso()函数，设置λ的值为1000\n",
    "lasso = Lasso(1000)\n",
    "lasso.fit(data.iloc[:, 0:13], data['y'])\n",
    "print('相关系数为：', np.round(lasso.coef_, 5))  # 输出结果，保留五位小数\n",
    "\n",
    "print('相关系数非零个数为：', np.sum(lasso.coef_ != 0))  # 计算相关系数非零的个数\n",
    "\n",
    "# 返回一个相关系数是否为零的布尔数组\n",
    "mask = lasso.coef_ != 0\n",
    "print('相关系数是否为零：', mask)\n",
    "\n",
    "data = data.iloc[:, 0:13]\n",
    "new_reg_data = data.iloc[:, mask]  # 返回相关系数非零的数据\n",
    "new_reg_data.to_csv('../tmp/new_reg_data.csv')  # 存储数据\n",
    "print('输出数据的维度为：', new_reg_data.shape)  # 查看输出数据的维度"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 模型训练"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# 自定义灰色预测函数 \n",
    "def GM11(x0):  # x0为矩阵形式\n",
    "    import numpy as np\n",
    "    x1 = x0.cumsum()  # 1-AGO序列\n",
    "    # 紧邻均值（MEAN）生成序列\n",
    "    z1 = (x1[:len(x1) - 1] + x1[1:]) / 2.0\n",
    "    z1 = z1.reshape((len(z1), 1))\n",
    "    B = np.append(-z1, np.ones_like(z1), axis = 1)\n",
    "    Yn = x0[1:].reshape((len(x0)-1, 1))\n",
    "    # 计算参数\n",
    "    [[a], [b]] = np.dot(np.dot(np.linalg.inv(np.dot(B.T, B)), B.T), Yn) \n",
    "    # 还原值\n",
    "    f = lambda k: (x0[0] - b / a) * np.exp(-a * (k - 1)) - (x0[0] - b / a) * np.exp(-a * (k - 2)) \n",
    "    delta = np.abs(x0 - np.array([f(i) for i in range(1, len(x0) + 1)]))\n",
    "    C = delta.std() / x0.std()\n",
    "    P = 1.0 * (np.abs(delta - delta.mean()) < 0.6745 * x0.std()).sum() / len(x0)\n",
    "    # 返回灰色预测函数、a、b、首项、方差比、小残差概率\n",
    "    return f, a, b, x0[0], C, P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "预测结果为：\n",
      "                x1       x3       x4       x5           x6       x7       x8  \\\n",
      "2014  8.14215e+06  7042.31  43611.8  35046.6  8.50552e+06   4600.4  18686.3   \n",
      "2015  8.46049e+06  8166.92  47792.2  38384.2  8.62714e+06  5214.78  21474.5   \n",
      "模型精度            好        好        好        好            好        好        好   \n",
      "\n",
      "          x13   y  \n",
      "2014  44506.5 NaN  \n",
      "2015  49945.9 NaN  \n",
      "模型精度        好 NaN  \n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "new_reg_data = pd.read_csv('../tmp/new_reg_data.csv')  # 读取经过特征选择后的数据\n",
    "data = pd.read_csv('../data/data.csv')  # 读取总的数据\n",
    "new_reg_data.index = range(1994, 2014)\n",
    "new_reg_data.loc[2014] = None\n",
    "new_reg_data.loc[2015] = None\n",
    "Accuracy = []  # 存放灰色预测模型精度\n",
    "l = ['x1', 'x3', 'x4', 'x5', 'x6', 'x7', 'x8', 'x13']\n",
    "for i in l:\n",
    "    f = GM11(new_reg_data.loc[range(1994, 2014), i].as_matrix())[0]\n",
    "    new_reg_data.loc[2014, i] = f(len(new_reg_data) - 1)  # 2014年预测结果\n",
    "    new_reg_data.loc[2015, i] = f(len(new_reg_data))  # 2015年预测结果\n",
    "    new_reg_data[i] = new_reg_data[i].round(2)  # 保留两位小数\n",
    "    C = GM11(new_reg_data.loc[range(1994, 2014), 'x1'].as_matrix())[4]\n",
    "    P = GM11(new_reg_data.loc[range(1994, 2014), 'x1'].as_matrix())[5]\n",
    "    if P>0.95 and C<0.35:\n",
    "        Accuracy.append('好')\n",
    "    elif 0.8<P<=0.95 and 0.35<=C<0.5:\n",
    "        Accuracy.append('合格')\n",
    "    elif 0.7<P<=0.8 and 0.5<=C<0.65:\n",
    "        Accuracy.append('勉强合格')\n",
    "    else :\n",
    "        Accuracy.append('不合格')\n",
    "\n",
    "new_reg_data = new_reg_data.iloc[:, 1:]\n",
    "new_reg_data.loc['模型精度', :] = Accuracy\n",
    "outputfile = '../tmp/new_reg_data_GM11.xls'  # 灰色预测后保存的路径\n",
    "# 提取财政收入列，合并至新数据框中\n",
    "y = list(data['y'].values)\n",
    "y.extend([np.nan, np.nan])\n",
    "new_reg_data.loc[range(1994, 2016),'y'] = y\n",
    "new_reg_data.to_excel(outputfile)  # 结果输出\n",
    "# 预测结果展示\n",
    "print('预测结果为：\\n', new_reg_data.loc[[2014, 2015, '模型精度'], :])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "真实值与预测值分别为：\n",
      "             y       y_pred\n",
      "1994    64.87    37.830236\n",
      "1995    99.75    84.400064\n",
      "1996    88.11    95.256053\n",
      "1997   106.07   106.863698\n",
      "1998   137.32   151.324295\n",
      "1999   188.14   188.310990\n",
      "2000   219.91   219.586734\n",
      "2001   271.91   230.376110\n",
      "2002   269.10   219.762444\n",
      "2003   300.55   300.550000\n",
      "2004   338.45   383.441026\n",
      "2005   408.86   463.106181\n",
      "2006   476.72   554.750685\n",
      "2007   838.99   691.243852\n",
      "2008   843.14   842.902195\n",
      "2009  1107.67  1087.404040\n",
      "2010  1399.16  1378.772344\n",
      "2011  1535.14  1536.534100\n",
      "2012  1579.68  1739.075209\n",
      "2013  2088.14  2085.603639\n",
      "2014      NaN  2187.909640\n",
      "2015      NaN  2539.105858\n",
      "预测图为： Axes(0.125,0.125;0.775x0.755)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEHCAYAAABP3uaxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl0VEXax/FvJWRjRwTCmqCioCiCUUFB4oKCOojACGMc\nlGUCuG8zLqgwbqMDbuOAGMVlxghuuDDqq6JEEQUJIigqIHswICISIAmEpN4/qkNCSEink+5Op3+f\nc/p0d/W9t6pzwpOibtVTxlqLiIjUbRHBboCIiPifgr2ISBhQsBcRCQMK9iIiYUDBXkQkDCjYi4iE\nAQV7EZEwoGAvIhIGFOxFRMJAvWA3oNiRRx5pExMTg90MEZGQsmTJkl+ttS0qO67WBPvExEQyMzOD\n3QwRkZBijNngzXEaxhERCQMK9iIiYUDBXkQkDNSaMfvyFBQUkJWVRX5+frCbEjSxsbG0a9eOqKio\nYDdFREJYrQ72WVlZNGrUiMTERIwxwW5OwFlr2b59O1lZWXTs2DHYzRGREFarh3Hy8/Np3rx5WAZ6\nAGMMzZs3D+v/2YjIYWRn0wWO8+bQWh3sgbAN9MXC/fuLyGFMnEh9aOjNobU+2FdFejokJkJEhHtO\nTw92i0RE/CAuDoyBZ57x+pQ6E+zT0yE1FTZsAGvdc2qqAr6I1EFr18If/lClUyq9QWuMaQLMAiKB\nPcAw4CdgreeQ66y13xpj/g5cCHxlrb3Gc+4hZb668Ub45puKP1+4EPbuPbgsNxdGj674j9/JJ8Pj\nj1enVSIiQdC6NSxfDoAF680p3vTsU4BHrbXnA1uA24GZ1tpkz+NbY8wpQG/gNOAXY8x55ZVV/Rt5\nr2ygr6zcWxMnTmTWrFkATJo06cBrEZGgmTvXDV+cfjor4QdvTjHWevVHwR1szOvA17g/AHuAb4Gx\nwPVAvrV2mjGmJzAA2Fm2zFo7scz1UoFUgA4dOpyyYcPBKR5++OEHunTp4lXbEhPddy8rIQHWr/f6\nKx5izZo13HTTTbzzzjv06tWLTz75hLi4ON8v6IOq/BxEpI7bv98NS+TmwvffY+Lillhrkyo7zesx\ne2NML6AZ8BFwnrX2NCAKN0zTANjsOfQ3oFUFZQex1qZZa5OstUktWlSatO2wHngA6tc/uKx+fVde\nHUcffTS7du0iIyODrl27BjzQi4gcJC0NVqyAKVMgNtbr07wK9saYI4AngVHAcmtttuejTKATsBso\njoINPdctr8xvUlLczyAhwd2kTkhw71NSqn/t4cOHM2rUKEaMGFH9i4mI+GrHDrjnHkhOhksvrdKp\nlQZgY0w08Bpwh7V2A/BfY0w3Y0wkMAhYBizBjc8DdAPWV1DmVykpbsimqMg910SgBxg6dCjGGHr3\n7l35wSIi/vL3v7uA//jjrldbBd70tkcDPYAJxpgMYAXwX+Ab4Etr7Vzgc6C7MeYJPDdwKygLOStW\nrGDAgAHceeedWuAkIsHz448wdSqMGQPdulX59CrdoD3shYyJAy4CvrbWrq2orCJJSUm27OYlujHp\n6OcgIlx4ISxYAKtXQ8uWB4qNMV7doK2xRGjW2jzg9crKRESkit5/3z2mTDko0FdFnVlBKyJSJxUU\nwM03Q6dOcN11Pl+mVqc4FhEJe9OmufH6d96B6GifL1P3evbZ2dC3L2zZEuyWiIhUz6+/wqRJ0K8f\nXHxxtS5V94L9fffB55/DvfcGuyUiItUzcSLs2gWPPVblqZZlhc4wTmWZ0ObPdxPsiz31lHtERECf\nPuWfU8syob3wwgsAXHXVVUFth4jUAt99B9Onw/jxcMIJ1b5c3enZn3aau0sd4flKERHu/emnB7dd\nIiJVZS3cdBM0aeIWUtWA0OnZe9MDHz/e5UiIjYV9+2DIEHdzoxomTpxIly5dGD58OJMmTaJz584M\nHz78oGMmTZrEokWLyM3NpUWLFsyaNYt69eqRnJzMJZdcwvPPP8/y5cux1pKamsqqVato0aIFr7zy\nCkVFRVx22WX8/vvvREVFcfnll1ervSJSB8yZ4zJb/utf0Lx5jVyy7vTsAbZuhXHjXHL7ceNq5Cbt\niBEjePnllwH44IMPuOSSS8o9rk+fPnz66ae0atWKt99+G4Ds7GyMMSz35J1+++23KSgo4NNPP6VD\nhw68++67zJ49m4SEBObNm0dCQkK12ysiIW7vXrjlFujSxcWxGhI6PXtvzJ5d8nrq1Bq5pLdZL085\n5RQATjrpJNZ7cio3adKE66+//sAxK1eu5MsvvyQ5OZndu3fTpUsXtm/fTjfP0uekpEoXwYlIXffk\nk/DTT/B//wdRUTV22brVs/cTb7JefvXVVwAsXbqUY445BoD69esTEVHyIz7uuOMYPnw4GRkZPP74\n4xx//PF06NCBFStWHDhXRMLY1q1uRuFFF8EFF9TopRXsveBN1svFixeTnJzM77//zsUVzIcdOHAg\nP//8M3379uWuu+4iISGBIUOGsGrVKpKTk1m1apW/voKIhIK77nKbkjzySI1fum4N4/jBihUrGDly\nZKVZL2+++WaSk5MPKsvIyDjofUREBM+UsyHuO++8UxNNFZFQtnQpzJjhppkfd1yNX17BvhInnHDC\ngSEa4JCA3qRJkwM3ZEVEfGKtC/LNm7vNSfyg1gd7a22tyiNftrfubzWVglpEarE33oDPPnMLQZs2\n9UsVtXrMPjY2lu3bt4dtwLPWsn37dmKrsM+kiISY/Hz461/hxBPdxiR+Uqt79u3atSMrK4tt27YF\nuylBExsbS7t27YLdDBHxh+xsOOMMt4/qxx9DPf+F5Fod7KOioujYsWOwmyEi4h+33+4CfWIinHOO\nX6uq1cFeRKROiotzwzfF1q93WS1jYyEvzy9V1uoxexGROmX/fnjxxZKtBYsnn9SvDykpsG6d36pW\nsBcR8bfCQkhPh+OPh6uuclMs+/cv6c3n50PjxhAf77cmKNiLiPhLURG88gp07QpXXOGGb958E5Ys\nca9rOHHj4WjMXkSkphUVucSMkybBihWuR//aazB4cMmeG35I3Hg46tmLiPggPd1NoomIcM/p6biV\nsG+9Bd27wx//6MboZ86E5cth6NCSQB8ECvYiIlWUng6pqbB3QzbzbF/yN2Tz5qg5/NbxFLj0Ujej\n5qWXXK9++HCIjAx2kzWMIyJSVRMmuOSUk7mXPsxnGd1otW8bGzYfxREvvOBm1vhxgZQvaldrRERC\nwA8b4oijZJ58K9wq/5b7f4YrrwxWsw5LwzgiIlV0Vtu1vM1AirN27SGOl0jhrHb+mydfXQr2IiJV\ndOPDrTkJt7d0PjHEspfcyMbc+JD/5slXl4K9iEgVpbT8iI6sZxGn05NFpDccR//uW0hJCXbLKqYx\nexGRqigoYO+4G8jiKL54IINv7owF/D9Pvroq7dkbY5oYY943xnxojHnTGBNtjJlhjPnSGHNXqeO8\nKhMRCWnTphGz9gf+FvkoV4wJnb0mvBnGSQEetdaeD2wBhgOR1tpewFHGmE7GmMHelPnrS4iIBMS2\nbdiJE5lXrx+Rlw48kM8sFFQ6jGOtnVbqbQvgCuBxz/sPgd5Ad+BVL8pWl762MSYVSAXo0KGDT19A\nRCRg7roLu2s31xQ9zr/G1p7tUr3h9Q1aY0wvoBmwCdjsKf4NaAU08LLsINbaNGttkrU2qUWLFj59\nARGRgFi6FJ55htmtr2XvUcf7e6+RGudVsDfGHAE8CYwCdgNxno8aeq7hbZmISOixFq6/nv1NmzNm\n8yTGjAlqmhufeHODNhp4DbjDWrsBWIIbkgHoBqyvQpmISOh55RX4/HPePvUB9tRrysiRwW5Q1Xkz\n9XI00AOYYIyZADwP/NkY0wYYAPQELDDfizIRkdCyZw/89a8Undyda74ezR/+4Nc9Rvym0p69tfYp\na20za22y5/EikAwsBM621u601uZ4U+avLyEi4jcPPwxZWXx66RNs/TWS1NRgN8g3xlpb+VEBkJSU\nZDMzM4PdDBGREuvXQ5cuMGgQ/X6dyerVsGZNrchYfIAxZom1Nqmy40LsFoOISADdeisYw/qr/8nc\nuTB6dO0K9FWhYC8iUp558+CNN+COO0h7vz0RETBqVLAb5TvlxhERKWv/frj+ekhMpOCGW3nuWLjo\nImjbNtgN852CvYhIWU8/Dd99B6+/zv8+jmPrVkL2xmwxBXsRkdK2b4e774azz4bBg0m70PXo+/cP\ndsOqR2P2IiKl3XMP7NwJTzzBho2GDz5wN2Zr2ZayVaZgLyJSbPlymD4dxo+HE09kxgxXPHp0cJtV\nExTsRUTA5b+54QZo2hTuvZf9+2HGDDd8UxeS8ob4f0xERGrIG29ARgZMnQpHHMH7c+Dnn93bukA9\nexGRvDy3gOrEEw9Mu0lLczlwLrooyG2rIerZi4hMngwbNriFVPXqkZUF770Ht98OUVHBblzNULAX\nkfCVnQ2DBrkbs0OHQnIyAM89B0VFMGZMcJtXkxTsRSR83XcffPWVS3gzZQoAhYXw7LPQrx907Bjk\n9tUgBXsRCT9xcZCfX/K+sBASEyE2lg9n57FpEzz6aNBa5xe6QSsi4WftWjjvvJL3cXGQkgLr1pGW\nBi1awMCBwWuePyjYi0j4WbPGTbMEiImBvXuhcWOybTxz5sDIkRAdHdQW1jgFexEJLwsWwIABEBsL\nV14JixbBuHGwZQvPP+9GdOrSjdliGrMXkfCxYIFbEtumjZtm2aaNK586laIiePYYl/+sU6fgNtMf\n1LMXkfDwxRcu0LdufXCg9/j4Y1i3LvRTGVdEwV5E6r4vvzxsoAe3YrZ5c7j00iC0LwAU7EWkbvvy\nS7jgAmjVygX6crab2roV3nrLDeHHxAShjQGgYC8iddfChS7Qt2xZYaAHePFFtxNhXbwxW0zBXkTq\npkWLSgJ9Rga0a1fuYdbCM89Anz7QpUtgmxhICvYiUvd89RWcfz4ceaTr0VcQ6MH9Hfjpp7p7Y7aY\ngr2I1C2LF5cE+owMaN++3MPS012GhHPOAWOgoCCgrQw4zbMXkbpj8WKXweyII1yP/jCBPjUVcnPd\ne2vh2mvdqtmUlAC2N4DUsxeRuiEzsyTQZ2SUu5fgrl3wySdwzTUlgb5Ybi5MmBCYpgaDgr2IhKT0\ndDi1XTafmr5c0fID9vbtB82auR59hw4UFcGPP8ILL8DYsXDSSW572XPPhZ07y7/mxo0B/QoBpWEc\nEQk5xcMwk3Pvow/zOX3bQLJNG2b1z2Dfiwl8+aWbjLNjhzu+SRM4/XQYPBh69oS//AWysg69bl3Y\nWLwixlob7DYAkJSUZDMzM4PdDBEJAfkRccTa/EPK84ilgcnj+OOhVy/36NkTOneGiFLjGGXH7AHq\n13eraENtzN4Ys8Ram1TZcerZi0jIOcqu4W0GksQSDC7Iv8EQ/soUduxwPfnDKQ7oEya4oZsOHeCB\nB0Iv0FeFV2P2xphWxpj5ntdtjTFZxpgMz6OFp3yGMeZLY8xdpc47pExEpDpWfZPL5IjbOJUlAOQR\nQzT7yKExMQnxlQb6YikpsH6922t2/fq6HejBi2BvjGkGvAg08BSdDjxgrU32PLYZYwYDkdbaXsBR\nxphO5ZX560uISN1XVATP3fkT+T168aeidL6nC9MYT08WMZ1xtI3YwgMPBLuVtZc3PftCYBiQ43nf\nExhjjPnaGPOgpywZeNXz+kOgdwVlBzHGpBpjMo0xmdu2bfPpC4hI3bdmDdx54hwG/yOJxHpZ7Hj5\nfZa+9D2TE6bxrenG5ISp7P7P7DrfO6+OSsfsrbU5AMaY4qL3gfuAXGCuMeYkXK9/s+fz34AeFZSV\nvXYakAbuBq2vX0JE6qaiIpg+tZCcmyfx0P772Z7Qgybz3sB0TCSFuj/0UpN8mWf/hbV2l7W2EFgK\ndAJ2A3Gezxt6rltemYiIV9avhyHJ2zn6+gu5ff/97Bk2iuY/LsB0TAxuw0KULwH4A2NMa2NMfeB8\n4DtgCSXDNN2A9RWUiYgclrVuCuSfj1/C45+fwnn1MrBPp9Fg1gy3b6z4xJepl38H5gH7gOnW2pXG\nmGxgvjGmDTAAN65vyykTEalQVpbLKd/mg+f4OOJqIlq3JPKtz+HUU4PdtJDndc/eWpvseZ5nre1s\nrT3JWvtvT1kO7obsQuBsa+3O8spquO0iEsKKs05GREBCgktpcMoJ+Vw2N5XnGE3U2X2ot+xrBfoa\nUmOLqqy1OyiZfVNhmYhI2RWsGzfC+2kbmRs7hBMLM+GOOzD33QeRkcFtaB2im6YiEnATJrhAH082\nGfRlKK/yNT1I2LsK3nwTHnxQgb6GKdiLSMAVZ5e8m3vpw3xeYRhbiOdUuxgGDQpu4+oo5cYRkYDL\nJY5YDk5k1pUVLDPdgLzgNKqOU89eRAJu8ri1vMcAildS7qE+syJTePfJdUFtV12mYC8iAbV9Ozz3\nZjN6mYUA5BNDHPn0SG7MkGvig9y6ukvBXkQCxlq4+mq47pe7aGZ3YC66iNhvFhFx9TiObbwl2M2r\n0zRmLyIBM3MmbHn1U24yj8K4cfDUU+6DqVOD27AwoGAvIgGxaRPcPn4nX0WPgA5Hw5QpwW5SWFGw\nFxG/KyqCkSPhwdwbaFWUhfnvAmjQoPITpcYo2IuI3/3739Do4ze5ghfhrrvcxrASUAr2IuJXP/wA\nj/5tC8uiU7Fde2DuuSfYTQpLmo0jIn5TUAB/vsLydNFfaByxG/Pf/0JUVLCbFZbUsxcRv7nvPuj+\n9bNcwP9g8uNw/PHBblLYUrAXEb9YtAhmPbCG5fVugr7nwnXXBbtJYU3BXkRq3J49cOUVhaRHXUlM\nTD14/nmXuF6CRsFeRGrc3/4Gg36azKksgBkvQfv2wW5S2FOwF5Ea9cEHsGDaN2RG3AND/giXXx7s\nJgkK9iJSg7Zvh3FX5fNhzBVENjvSpUMwJtjNEhTsRaSGlE5y1qloBTz3HjRvHuxmiYfumIhIjTiQ\n5Mx6kpwNGBDsJkkp6tmLSLVlZcHtV+ewKOZKaK8kZ7WRevYi4pP0dEhMdDMqO3WC+3bdQHzBJrdK\nVknOah0FexGpsvR0SE2FDRvcWH3//De5sugFvvvDnUpyVktpGEdEqmzCBMjNhXiymc1gOrGKJfRg\n2NJ7+CnYjZNyKdiLSJVt3Oie7+ZeerKQQiI5i/+ydpOSnNVWCvYiUmW5No5Y8g+8r0ch33MC+cQC\necFrmFRIY/YiUmV3DFvLJyRjPe/3UJ9ZkSm8++S6oLZLKqZgLyJVsmQJfDx7J71ZAEAescSRT4/k\nxgy5Jj7IrZOKKNiLiNd27IDRg3fwlh1IVLTBXHEFcd8sJOLqcRzbeEuwmyeHoTF7EfGKtTD6yv1M\n2XQZiZHrMZ/MgzPPdB9OnRrcxkmlvOrZG2NaGWPme15HGWPmGGMWGGNGVaVMRELXI49A3zm3cJ6d\nS0Ta0yWBXkJCpcHeGNMMeBEoXhJ3HbDEWnsmMNQY06gKZSISgubPh9V/e4Yb+Bf2xptg5MhgN0mq\nyJuefSEwDMjxvE8GXvW8/gxIqkLZQYwxqcaYTGNM5rZt26reehHxu19+gccu/Yx/26spOK8/ZvI/\ng90k8UGlwd5am2Ot3VmqqAGw2fP6N6BVFcrKXjvNWptkrU1q0aKFb99ARPymsBBuunQ9T28fQmHi\n0US9NhPq6VZfKPJlNs5uIM7zuqHnGt6WiUgIeWjCLm77YiCN6+8n9sM50LRpsJskPvIlAC8Bente\ndwPWV6FMRELEB+8XccLDf+YE8z3Rb77qUltKyPLl/2MvAu8ZY/oAxwOLcMM13pSJSAjYtAlWDLmb\nm3mbfZP/ReT5/YLdJKkmr3v21tpkz/MGoB+wADjPWlvobVlNN15Eal5BATx77kxuznuQ3//4F6Jv\nvjbYTZIaYKy1lR8VAElJSTYzMzPYzRAJe49dvphxM89iV5fTaPnNRxAdHewmyWEYY5ZYaw+Z7ViW\nbpqKyAHvPfszl80cRG6jeFp++roCfR2iOVQiAsCa7/JoNXYQzSJ2Uu+TL0DToesU9exFhLxcy8qz\nxtC9KJPd09OJTjop2E2SGqZgLxLCSm/6nZjo3lf1/FPbZZPd4Ggu3PEy/9f7flr+5RJ/NFWCTMFe\nJESV3fR7wwb33tuAX3z+3zePpiPrWMUx/HHJHVX+gyGhQbNxREJUYqIL8GU1aACXXOKmUBYUwL59\nJa9Lv//imzjiSm0tWCzfxBJbpK0FQ4Vm44jUccWbfpe1Zw8sWgTffgtr1sDWrbBrl8tzExsLzZtD\nQtv9pJEKQBHGnUd9XiKFjlZbC9ZFmo0jEqLatoWsrEPLExLgp58Oc+LPP8Of/gR8xvd05jhWkUcM\nseSTQ2NiErS1YF2knr1IiOra9dCy+vXhgQcOc9KHH8LJJ0NmJl+M+w+rI7ownXH0ZCHTGUfbiC2H\nP19CloK9SAhaswY+/hj69nU9eWPcc1oapKSUc0JhIdx9N/TvDy1bQmYmZzz1Z3b/ZzaTE6byrenG\n5ISp7P7P7PLPl5CnYRyREHT77RAVBTNnQuvWlRycnQ2XXw4ZGTBqFDz5pPsvAO4Pg4J7eFCwFwkx\nCxbA66/DpEleBPq5c100370bXnwRRowIRBOlFtIwjkgIsRZuuQXatIFbbz3MgYWFMHEinH8+HHkk\nLF6sQB/m1LMXCSGvvuqmVT73nJtPX64tW9ywzbx5cOWVMHXqYQ6WcKFgLxIi9u51Y/UnnXSYTvrH\nH7thm5wceP55uOqqQDZRajEFe5EQ8eSTsH49fPQRREaW+iA7G4YNg9NOg0cfhc6dXdA/4YRgNVVq\nIQV7kRCwfTvcfz8MGADnnVfmwzvugPnz3WPECJg2TcM2cggFe5EQcO+9LuXB5MmlCuPiIL9Mbpv/\n/McN7Ocpt40cTLNxRGq51atdZ33MmFIjMzt3wqBB7rVxuW2oX9+N169Tbhs5lIK9SC13220ugdm9\n93oK5s6FE090Pfju3V2wj411vfzGjSFeuW3kUAr2IrXYZ5/Bm2+6WTitGuyGq6+Gfv1cL/6LL1ye\n43HjYOFC97xlS7CbLLWU8tmL1FJFRdCzp0tS+dNznxE7fqQbornpJne3Ni4u2E2UWkD57EVC3KxZ\n8O3iPD7sejOx/ZNd4aefwiOPKNBLlSnYi9RC+fnwys2L+D6mO8d/8BiMHw/LlkGfPsFumoQoTb0U\nqW327mXpgEnM3vpPClq2hfSPyplcL1I16tmL1CZff83+7kn0yniIj9tfReyqbxXopUYo2IsEW3Y2\nnHWWS2N5+uns2bCdgRH/o/0HM6BJk2C3TuoIBXuRYLvxRpfq4JFHyBkwjGPyv6Pd2Ivo0iXYDZO6\nRGP2IsESEwP79h1U1HhOOpt4g5xJSncgNUs9e5FAstZt+t2vnwv0kZFQz/W5CmPq8xIpPHPnOlq2\nDHI7pc5RsBcJhH37XJKybt3gggtgxQp46CGXpbKoCBsbi9mbj23YmDF3Kd2B1LwqB3tjTD1jzEZj\nTIbncaIx5u/GmMXGmKmljjukTCTs7NzpUlUedZTbNaqoyG0qsm6dS3rz++8wbhzv3rWQpxhH385b\ntF5K/MKXMfuTgJnW2tsAjDGnAL2B04B7jDHnATvKlllr59ZQm0Vqv02b4IknIC3N5SY+91x49lnX\nqy/OUgmkD5nNHXe4w6Ojp9LkBrgiiM2WusuXYN8TuNgYczbwLbASeMNaa40xHwADgJ3llB0S7I0x\nqUAqQIcOHXz8CiJBlp0Nw4fDK6+4RGRTprjX1rodpG65BXr0OOS09HRITYXcXPd+3z4YO9b9LUhJ\nCfB3kDrPlzH7xcB51trTgCggDtjs+ew3oBXQoJyyQ1hr06y1SdbapBYtWvjQFJFa4N573dTJpCSX\ncvjtt+G662DtWhfRywn0RUXub0BxoC+WmwsTJgSo3RJWfOnZL7fW7vW8zqQk4AM0xP0B2V1OmUjd\nEhvrdgEvttnTvykocHvBllFUBIsWuU7/66/D1q3lX3bjRj+0VcKeL0H4v8aYbsaYSGAQrhff2/NZ\nN2A9sKScMpHQZ23JXq/Fisfg4+Lc+Mv69QcdvmiR68UnJsIZZ8D06XDqqXDkkeVXoRFN8Qdfevb3\nAi8DBngHuB+Yb4x5AujveWwA/lGmTCR0/forvPiiu8n644/QqBGMGsWGzG20XzybvcQQk7eXn7Y0\nplOreJZkuo2kXn0VNmyA6Gh3b/bBB2HgQLehVNkxe3B7kjzwQPC+ptRdNbJ5iTEmDrgI+Npau7ai\nssPR5iVS6xQVwbx58MwzMHu2G57p1Qv+8he47DLS32pAwxGDySpqTRqppJJGG5NNavPZ/PqrWyt1\n/vnuHu3AgdC06aFVpKe7MfqNG12P/oEHdHNWqsbbzUu0U5VI6dk08fHu/QsvuF782rXQrJkbthkz\nBrp2PXBaYqLrtZcVG+s2CB80yJ0q4k/eBnvlxhG57z74/HMYPdqNt8yZA4WF0Levm2kzZIiL4B5r\n17pdpMoL9ODu2Y4cGaC2i3hJwV7CV1yc2xKq2HvvuefISFi5Eo499sBHW7a48feXX3Y3XMHlMSs9\nGaeYbrBKbaQpkRJ+iopg7lw4++yDy6Oj3XBOVhYceyy//w7PPedylrVtCzfc4P42PPywm3AzY4a7\noVqabrBKbaWevYSP335zY/HTp8Pq1dC8Ob916EbTjcvdbJp9+1i5tRnffhbPzJmuo79vHxx9tLuJ\n+qc/cVCO+YQE96wbrBIKdINW6jZr4auv4Kmn3A3Y/Hw3o+bqq5lZMJT6Yy4/aDZNPNkMZTatW7tO\n/p/+5BbGlkpnI1Kr6AathLc9e2DmTDctZulSaNAArrySvaPGs6JeN5YuhZtugl1Fsw+cci0uQWvL\nli4xWWRksBovUvMU7CWkpafD47dlM2XzcG5t9wr3XLeDP2Q95XLH79zJ7o5dWTR4Kq9EXcHCLxrz\nwwzYv//w19y2TYFe6h4FewlZxStQH8mdRB/m81ZWD9relk2BieLduKFM4WoWrDsT1hnatoWTT3aL\nm04+2T3OPbf8PDSaTSN1kWbjSLWkp7vFRRER7jk93X/n5+fDd9+5JGL33w9Dr4hlT65hHGlEYGlL\nNgCFNoLinwZGAAAOsElEQVTXBr3MJf/szYcfGn75xU2w+d//POcNhWOOcakLNJtGwoV69uKzsrld\nNmxw78G7GSkVnb9rF5xwgpvq/uOPJY916yCiqICBvEMqacSwlyKgiEjqUcge4niTwfyVKWR78Uen\nuI2aTSPhQLNxxGft2pVk9S0tLg4uvtjlhineT7v4ufTrZ5+FnJzD1xEb69Y29Wm3jst2PsOp3z1H\n3M6tFLVrz+O7xnDkztWk8DL7iCaafTzNWP6ZMK104kmROk2zcaTGZWe7vGCffOIexYE+nmxmMZxh\nvMJW4snLg2+/dRkH9u93j+LXpcv27Km4rvffh+OOKqDD8jlEPvM0vP+Rm/940UUwdiwR/fvTalYk\nDUcMZnrRuANTJ9tFZGsYRqQ81tpa8TjllFOsBN5LL1mb1PZnm8FZNqldtn3ppZLPtm+39o03rL3m\nGmu7dLHWTVq3tmlTawcNsrZZM/d+KuPtfiLsVMZbsDYhwbu6ExLc+fG4+luRbcHaM9uus/bOO62N\nj3cHtGtn7aRJ1m7cWG77ExKsNcY9l26/SDgAMq0XMVbDOGGseMx8cu7VjOVpnmYsN8dM49xzXS/+\nm29ceG9Q39LvjD30Oz2Hs7rn0KXtLiL35FB4wQAiCwsOuW5hZDSRCxe4nO+NG7vnBg0OWZlUtv4P\n6Ue9CDjPfogxBi680G3K2r+/G/cRkUMoxbEc1vY1v9PwmFbEsO+Qz4owbKzfhWb1cmhQmENk7i5M\ndX9PIiKgYcOS4N+oEWRmujw1ZdWr51JLtm9fvTpFwoDG7MNE2UVFNz4Uf2A2SeG+QjZ/uZGtn/7I\n7q9XErHyRxr//CNtdq2kld1y4BoWt+1YEYbtHMF3dOXsAUe6wFz8KO6ll30/ebLL9xsV5RLJDBkC\n117r7rzu2nX45+7d4aef3GtrXRrJwYPd/q3x8UH5eYrUVQr2Iax4GOSx3Hvow3wezLqCr0acwdzr\nf6RD7o+0z19NB/IpXiP0mzmCzY06s7rThazsdBzPft6Z/jtnMoxXD8xmeY3L3GyW171sRF4ejBvn\nGpKW5sZ/+vb1/kuMH+/Oi4lxfyyaNlWgF/EDDeOEEmth0yYKMpex+b1ltJ8xkUgOHQYpwpDZ8iJy\nO3Qm8vjjaHJ6Z9qd15kjjj14h+v0dA7ZVq9dRDa7/zM7cHPNBw+G1q0P/mMxe3bl54kIoDH70FJ2\nWzxwPeYVK2D5cgq/XsbuBcuIXrmcuLwdB05bRwKRFBLPVqIpII9Y3mQQt/AY2da73rH2QBUJbRqz\nDxHp6RA19k6G7pnPNwmDOKJHRxJ+XwarV2EKCwHIowHfcyLLuIxf23aj4RndOObSExn/t0bcmTWe\nVNLII5Zo9vE7zYhJ8H4YJCVFwV0kHCjYB1pBASxbBgsXUnT9DaTYkmGYHvsWwcJFFBLBg9zJMrqR\n16kbx1xwNH3PjuCPfaF585JL7SyChiO2alGRiFTOm8n4gXiE6qKqwy1KstZau3mzW5l066123+m9\n7f7o2AOrk7bQ0m6gnd1LlLVg9xBn/0uKPap+tn3tNWt/+cW7+rWoSCR84eWiqpDPelndrItAyQyS\nLVsqP7ZM3ampMHLzffTmc0ZnTeSZ0QuZd8njbO49jJ3NEtzmpUOGsHfKv8hctJ9/7RvHZbzCMVEb\niGcL73IxkRSSRywx7CWHxqzLi2foUGjRovI2pKS4/VCLityzhmREpDwhfYO2ONg2zi3JzbKrfjxp\nad4HvfR0KBp3NZfvfpqXG44lYvo0d+7+/W7+986dBx77fs0hZ9NOdm/eSd6WnXR66R7qUVjudTfQ\ngS/pxYqGPfn1mJ5EnNKdo4+PoXNn6NzZ7V969NHwyIbBbOHgbfFuSZitRF4i4pWQm41jTJJNSMg8\n7GyQ4t7rd9+5RFsPPujS406lZLn/TTxGE3I4rnUO8Q120SouhyOj3eOIejk0jdxFE5NDQ5vDiZ9O\nJcIeGqyLFxlVZh/1KCKCaAqIwFJAPRZwBjfwBE9/eTLHHQfNmlV8ftkUv+DyqVflj5WIhLeQDPaQ\neSDYnXOOC+rFj2+/dTMRiwNjU3awlVZEc2hulsoUEkEOjdlDfRqyh0bsIpIi9hPJOhL5iH78HtcG\nGjcmomkToo5sQmzLxsTFN6Fh2yY0ad+YIzo24dI/xTJhs5sN42uKXU19FJHqCNlgD278vXTKlJYt\nLOcfs5b+jRbQI38BCVkLqL92BQBFgMUQ6elZf8cJfNh0GLc90bbC5f75EfXZmWNo3Rqm2kOD9bVm\nWrkpW8qqFYuSRCSshVywP840sjtZzVbiqUcB6bcspXvuAtpvWkBs5oKSm6dNmkCvXnDmmczNO5NN\nD6czovD5A8F6RuRYGrw4zatgm5hY/TFz9cxFJJhCLtgnGWOn04PdNOJ08xVxNs99kJgIZ54JvXu7\n5xNOcF1/j42nDibjx9Y8ujuVmxumkdw5mw6LvVturzFzEQl1IRnsi+fiFEVEEjFrpgvubdr4tV71\nzEUklNWaYG+MmQEcD7xrrb2/ouOSjLHzTRxbeg2m4xtTlPlQRMQL3gZ7vy6qMsYMBiKttb2Ao4wx\nnQ5zMHFmLx27NVagFxGpYf7OjZMMvOp5/SHQG1hd7pFdukByslvNKiIiNcrfwb4BsNnz+jegR+kP\njTGpQCpAhw4dYOpUPzdHRCQ8+Ts3zm4gzvO6Ydn6rLVp1toka21SC28SwYiIiE/8HeyX4IZuALoB\n6/1cn4iIlMPfwzhvAfONMW2AAUBPP9cnIiLl8GvP3lqbg7tJuxA421q705/1iYhI+fy+U5W1dgcl\nM3JERCQIas0KWmPMLmBlEJtwJPCr6lf9qj+s6q4L9SdYayud4VKb9qBd6c0qMH8xxmSqftWv+sOr\n7nCqP+S3JRQRkcop2IuIhIHaFOzTVL/qV/1hWX84f/eA1V9rbtCKiIj/1KaevYiI+ImCvYhIGAhY\nsDfGtDLGzPe87mGMmWuMWWCMuaWislLndjXGfBTE+ucYY04OdP3GmKOMMR8bY74xxtxdnfpFJMxZ\na/3+AJoB/wd87Xm/AGgPGOALoGN5ZZ5jDS4XfkaQ6k8BHg/G9wceBc70nPM50MKHupsA73t+hm8C\n0cAM4EvgrlLHeVUWqPrLOy+Q9ZcqbwUsDfT3L1U+DfhDEH7+zYD3gEzg6QDU3wqYX+p9FDDH8+9i\nVBDq7wBkAJ/gbqCaQNZfqrwr8JGv37/0I1A9+0JgGJDjeX+EtXaTdd9mO9C4gjKAkcC8YNRvjDkC\neATYYYw5O9D1e55PMsa0AmKA332oOwV41Fp7PrAFGE6Z3cPK21GsSruM+aH+cs7rH+D6i02hJE13\nQOs3xvQB4q21c4JQ/5+BdOsW+zQyxvi66Meb+psBL+L2vyh2HbDEWnsmMNQY0yjA9Y8Fxltrz8F1\nwk4McP0YYwyuwxflY90HCcgKWusSouHaDsACY8y1uA1NEoHl5ZUZY5oDVwAXeB4BrR+4F3gNeBr4\nhzGmkbX2nQDWXw+4HmiH62Hs96HuaaXetsD9PB/3vC/ePaw7h+4oVl5Z+buM+aH+cs77pap1V6d+\nYLUx5hxgD+4fqU98rd8Ysx54BnjPGHOJtfbtQNaP62h0NcY0xQW7TX6s/w1cZ6j0d0wGbve8/gxI\nwodOn6/1W2snlDqvOT6mM6jG94eSjq7Psa+0YN2gHQv8CFwLPOzpzZZX9hBwh7W2IEj1dwemWmu3\n4P4xJAe4/tuBqzy/eHFAP18rNMb0wv3XfBMH7x7WikN3FKuozGc+1H/QedbahYGs3xgTDdxNScCp\nFh++/wjge+CfwGnGmOsCXP/nQAKus/GDp9wv9Vtrc+yhGXED9vtXQf3F5w0DVlhrfw5k/aU6ulOq\nU29pQQn21tpCSpKepVdUBvQFHjbGZAAnG2PuD3D9PwFHeV4nARsCXH9HoL0xJha3paNPiyI8w1FP\nAqMof/cwb8t84mP9Zc/zmY/13w5Ms9b6MnRWE/V3B9I8HY2XAJ+HEX2sfyIwzlp7L64TMtKP9Zcn\nkL9/FZ13FHArcKOvdVej/hrv6AZz6uX9wG2eHmy5ZdbaY621ydbaZOAba+1dgawf16u61hizADgL\neC7A9U/E3STahusRfFLVSjw91NdwvzgbKH/3MG/LqszX+ss5zyfV+P7nAdeU6mg8G+D6a6SjUY36\nmwEnGmMigdPxvaPhTf3lCeTvX3nnNQNm4m4O+7wPRzW+f813dL25i6tH6D6A8cAO3B+NDOBKYBnu\nxs8PuNkCjb0pC3D9Zc8bFsj6y1wjIwg//0a4IPEZbuZG2wDXfxqwAtcT/Qho6K/6y/s544aQVgBP\nAItxNzUDWf/DQHap8/oGsv6a+v0r/VC6hDDk6bX0Az6zbpjA67JA1u8vqr/21V/BcW1wveAPbA3u\nchfI71qb6lewFxEJA0qXICISBhTsRUTCgIK9CGCMiTHGHB3sdoj4i4K9iHM5bi50hYwx/zDGJBlj\nIowxfY0xRxtjxgSofSLVohu0EvaMMS2Ab4A1uPxFrXHT7iKAOGvt2Z6FbR/hVlH3BK4CrgbetS7v\niUitFpDcOCK1lTEmCrdi+Ulr7UOesoXW2ovLHDoW+NRaW2iMGY9LaVFgjFlpjDnDWvtFgJsuUiXq\n2UtYM8Yk4paxNwaO8RSfiVvIVA94B5eg6gtcatpPgZHW2hGe85sDs4GLrbW7Atl2kapQsBcBjDGf\nAOdba/d7evY9S312OXAkcDQudcBm4BzcSsj2uJWe71prnw58y0W8o2EcCWue3C9QQe4XY0wE8ArQ\nCzgCF+T3Ae9Yay80xtwBZFprq7WTmoi/KdhLuEsFLgHygLc8ew4cZ4z5n+fzSNwOQ98DWGtzjTGn\nAt95Po8BcgPaYhEfKNhLWLPWPgU8VbrMGLOo7A1aY8xZQIQni+EkXOpbcBtS1FjeFhF/UbAXOVTD\ncspicPuHPgK8bK39wRjzAi4V8KoAtk3EJ7pBK1IFxhhj9Y9GQpCCvYhIGFC6BBGRMKBgLyISBhTs\nRUTCgIK9iEgY+H9AfeSWYCmLjwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x203551f2940>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.svm import LinearSVR\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "data = pd.read_excel('../tmp/new_reg_data_GM11.xls')  # 读取数据\n",
    "data = data.drop('Unnamed: 0', axis=1, errors='ignore')  # 如果该列存在则删除，不存在则跳过\n",
    "data = data.drop(labels='模型精度', axis=0, errors='ignore')  # 删除行\n",
    "feature = ['x1', 'x3', 'x4', 'x5', 'x6', 'x7', 'x8', 'x13']  # 特征所在列\n",
    "data_train = data.loc[range(1994, 2014)].copy()  # 取2014年前的数据建模\n",
    "data_mean = data_train.mean()\n",
    "data_std = data_train.std()\n",
    "data_train = (data_train - data_mean) / data_std  # 数据标准化\n",
    "x_train = data_train[feature].as_matrix()  # 特征数据\n",
    "y_train = data_train['y'].as_matrix()  # 标签数据\n",
    "linearsvr = LinearSVR(random_state=123)  # 调用LinearSVR()函数\n",
    "linearsvr.fit(x_train, y_train)\n",
    "\n",
    "# 预测2014年和2015年财政收入，并还原结果。\n",
    "x = ((data[feature] - data_mean[feature]) / data_std[feature]).as_matrix()\n",
    "data[u'y_pred'] = linearsvr.predict(x) * data_std['y'] + data_mean['y']\n",
    "outputfile = '../tmp/new_reg_data_GM11_revenue.xls'\n",
    "data.to_excel(outputfile)\n",
    "print('真实值与预测值分别为：\\n', data[['y', 'y_pred']])\n",
    "\n",
    "print('预测图为：', data[['y', 'y_pred']].plot(style = ['b-o', 'r-*']))  # 画出预测结果图\n",
    "plt.rcParams['font.sans-serif']=['SimHei']    # 解决中文乱码\n",
    "plt.xlabel('年份')\n",
    "plt.xticks(range(1994,2015,2))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 性能度量"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "平均绝对误差： [ 34.24652892] \n",
      " 中值绝对误差： 17.8079477946 \n",
      " 可解释方差： [ 0.99086334] \n",
      " R方值: [ 0.99085445]\n"
     ]
    }
   ],
   "source": [
    "from sklearn.metrics import mean_absolute_error  # 平均绝对误差 \n",
    "from sklearn.metrics import median_absolute_error  # 中值绝对误差\n",
    "from sklearn.metrics import explained_variance_score  # 可解释方差\n",
    "from sklearn.metrics import r2_score  # R方值\n",
    "import pandas as pd\n",
    "\n",
    "data = pd.read_excel('../tmp/new_reg_data_GM11_revenue.xls')  # 读取数据\n",
    "data = data.drop('Unnamed: 0', axis=1, errors='ignore')  # 如果该列存在则删除，不存在则跳过\n",
    "mean_ab_error = mean_absolute_error(data.loc[range(1994, 2014), 'y'], data.loc[range(1994,2014), 'y_pred'], multioutput = 'raw_values')\n",
    "median_ab_error = median_absolute_error(data.loc[range(1994, 2014), 'y'], data.loc[range(1994, 2014), 'y_pred'])\n",
    "explain_var_score = explained_variance_score(data.loc[range(1994, 2014), 'y'], data.loc[range(1994, 2014), 'y_pred'], multioutput = 'raw_values')\n",
    "r2 = r2_score(data.loc[range(1994, 2014), 'y'], data.loc[range(1994, 2014), 'y_pred'], multioutput = 'raw_values')\n",
    "\n",
    "print('平均绝对误差：', mean_ab_error, '\\n', \n",
    "      '中值绝对误差：', median_ab_error, '\\n', \n",
    "      '可解释方差：', explain_var_score, '\\n', \n",
    "      'R方值:', r2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
